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The bold diagrammatic Monte Carlo (BDMC) method performs an unbiased sampling of Feyn- 
man's diagrammatic series using skeleton diagrams. For lattice models the efficiency of BDMC can 
be dramatically improved by incorporating dynamic mean-field theory solutions into renormalized 
propagators. From the DMFT perspective, combining it with BDCM leads to an unbiased method 
with well-defined accuracy. We illustrate the power of this approach by computing the single-particle 
propagator (and thus the density of states) in the non-perturbative regime of the Anderson local- 
ization problem, where a gain of the order of 10"* is achieved with respect to conventional BDMC 
in terms of convergence to the exact answer. 

PACS numbers: 02.70.Ss, 05.10.Ln 



A skeleton diagrammatic series is nothing but Feyn- 
man's diagrammatic expansion in terms of 'dressed', 
or 'bold-line', propagators, interaction lines, and ver- 
tices, which account for the summation of certain 
subclasses of diagrams. Its power lies in the fact 
that, even when truncated to the lowest orders, it 
often captures the basic physics of strongly correlated 
systems and yields quantitatively accurate answers. 
Among its numerous successful examples we mention 
screening effects, self-consistent Hartree-Fock schemes, 
the GW-approximation for simple metals, Bogoliubov 
and Gor'kov-Nambu equations, etc. Often, as, e.g., in 
case of Kohn-Sham orbitals in density functional theory, 
the diagrammatic structure is hidden in a set of integral 
equations, whose implementation has been improved to 
perfection. Physically, the lowest-order skeleton graphs 
embody the idea of incorporating some 'mean-field' 
theory self-consistently. 

The notorious shortcoming of self-consistent treat- 
ments based on the lowest-order diagrams is lack of ac- 
curacy and control: the error due to truncation can be 
established only by reliably calculating contributions of 
higher-order diagrams, which in the typical case of opti- 
mized codes solving a set of self-consistent integral equa- 
tions is nearly impossible (in the absence of small pa- 
rameters order of magnitude estimates are essentially 
meaningless) . The recently developed bold diagrammatic 
Monte Carlo (BDMC) method pi allows one to sample 
skeleton Feynman's expansions far beyond the mean-field 
level. Given that even the diagrammatic Monte Carlo 
method based on bare propagators can produce very 
accurate results for correlated systems (say, for the re- 
pulsive fermionic Hubbard model [5]), BDMC emerges 
as a powerful generic field-theoretical method. It has 
been successfully applied to the fermi-polaron problem 
[1], and, very recently, to the problem of equation of state 
in a system of resonant fermions [1] . The above examples 



deal with continuous-space problems, but it is natural to 
expect that working with the skeleton series will bring 
significant advantages to lattice models as well. 




FIG. 1: Schematic representation of the BDMC-I-DMFT pro- 
tocol (see text for the details). 

In this Letter, we show that, in addition to simply go- 
ing from bare to skeleton expansion, a dramatic increase 
in performance can be reached by employing an exact 
series re-summation procedure which accounts for the 
summation of all local contributions to the self-energy. 
This approach amounts to embedding the dynamic 
mean-field theory (DMFT) [5] solution into an exact 
diagrammatic method and avoids, in particular, any 
double counting or other uncontrollable errors. The gain 
in efficiency comes from two related observations: an 
impressive success of DMFT applications [51 [5] , and the 
fact that summation of local contributions can be done 
separately by a variety of highly efficient methods. The 
BDMC-I-DMFT approach thus involves two distinct but 
cross-linked numerical processes: (i) a problem-specific 
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solver of the DMFT-type problem (to be referred to 
as 'impurity solver', in accordance with terminology 
accepted in literature), and (ii) a generic BDMC scheme 
simulating skeleton diagrams which cannot be reduced to 
the purely local ones. The protocol is illustrated in Fig.jl] 

Below we start with the precise formulation of the 
combined scheme, and then proceed with its implemen- 
tation for finding a disorder-averaged single-particle 
propagator (and thus the density of states) in the 
non-perturbative regime of the Anderson localization 
problem which is well suited for illustrating the idea 
because the efficiency gained by incorporating DMFT 
solutions within the BDMC is about 10^. In general, the 
gain will be problem and parameter specific, and will 
also depend on the efficiency of the impurity solver. [We 
stress that our goal in this article is to explain the new 
method and illustrate its implementation, not to solve 
the localization problem in its full complexity.] 

Formalism. - The protocol of reformulating skeleton 
series to account for all local contributions to self-energy 
is conceptually straightforward. The Dyson equation re- 
lates the Green's function, G, to the self-energy, E (for 
clarity, we suppress below the frequency variable): 



G(p) 



1 



Goi(p)-S(p)' 
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with Go standing for the non-perturbed Green's function. 
The local propagator Gioc is defined by integrating over 
the Brillouin zone 'BZ' 
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We now separate contributions to the self-energy into two 
parts 



S(P) 



S'(P) 



(3) 



where Sjoc is given by irreducible skeleton diagrams 
which involve exclusively Gioc propagators. In other 
words, this local propagator has only purely momentum 
independent building blocks, while all the rest is put in 
S'. 



Numerically, one calculates the self-energy using 
current knowledge of the Green's function and then uses 
it to permanently improve the knowledge of G within the 
self-consistent process. This involves two steps. First, 
the current knowledge of Gioc serves as an input for the 
calculation of Eioc = Sioc[Gioc] achieved by the impurity 
solver, and G(p) is used for the BDMC simulation of the 
remaining skeleton graphs. Second, self-energies Sioc 
and E' are combined into the total self-energy, Eq. ( 3 1 , 
which is then used to find the updated G by Eq. ( 1 1 . 



This is illustrated in Fig. [T] 

Technically, the crucial advantage of separating local 
contributions to the self-energy is that the corresponding 
momentum independent problem admits a variety of 
techniques for solving it very efficiently 7\. Treating the 
local physics non-perturbatively is very appealing from 
the physical viewpoint. In typical problems such as the 
Hubbard model, the diagrammatic technique expands 
around the non-interacting limit which is dominated 
by large hopping processes. The competing phase 
with large on-site interactions tends on the contrary 
to localize the particles. Hence, building diagrams on 
top of the solution capturing essential physics of the 
competing phase may be better suited for describing the 
difficult intermediate regime as well. Local physics is 
also dominant at high temperatures which can easily be 
understood in terms of Feynman's path integrals. 

From Eqs. ([l])-([2]) it is explicitly seen that 
BDMC-I-DMFT process builds an exact solution of 
the problem on top of the DMFT answer, which is 
crucial not only for improving the quality of the final 
result but also for reliable estimates of corrections to 
mean-field results. 

One of the solvers for obtaining Sioc in terms of Gioc 
widely used in the standard DMFT approach is based 
on an implicit formulation of the problem in terms 
of the single-site (or impurity) effective action with a 
certain auxiliary (to be determined) 'bare" propagator 
go- The advantage of this formulation is in the flexibility 
of designing efficient tools (impurity solvers) for 
obtaining the Gioc [50] relation; the local self-energy 
readily follows from Eioc = ^ ^\oc- Iterations 
leading to the self-consistent solution consist of plugging 
the thus obtained self-energy in Eq. ^ to redefine the 
auxiliary propagator by ^q"^ = GjqJ[S] + Sioc- Solvers 
based on the effective action approach play a crucial part 
when the diagrammatic expansion of Sioc[Gioc] cannot 
be used because of technical or convergence problems. 

Illustration. - We illustrate the introduced concepts 
for Anderson's model of particle localization on a disor- 
dered three-dimensional cubic lattice. We consider delta- 
correlated gaussian disorder in the chemical potential, for 
which the standard diagrammatic technique can be for- 
mulated [5j. The Hamiltonian, in standard lattice nota- 
tion, reads 



(i-j) i 



(4) 



The random on-site potential is distributed with the 
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gaussian probability density 



P(e) 



(5) 



the dispersion V characterizing the strength of the 
disorder. We choose J = 1 as our unit. We work in 
the real-time representation where the Green function is 
defined as G(r,i'-i) = -i ( rc(0, t) c^(r, t') )■ We took 
a lattice of size 12 x 12 x 12. Just like in conventional 
DMFT, larger lattices pose no problem at all; in fact, 
larger lattices would suppress revivals and make hence 
the simulations easier. The (local) density of states is 
given by the imaginary part of its Fourier transform 
for r = 0, DOS(w) = -7r-ilmG'(r = 0,u;), which can 
be compared with the exact diagonalizaton results of 
Refs. [illl]. 

Evaluating the sum of all skeleton diagrams involviong 
local propagators only (i.e., the DMFT part (TU]) simpli- 
fies for Anderson's localization since disorder lines have 
no time dependence. For a single-site problem, one does 
not even need to expand the gaussian exponential into 
the diagrammatic series, because averaging the Green's 
function — in the frequency representation, the former 
is immediately found to be equal to l/[l/^o('^) ~ — 
over the disorder amounts to performing a simple one- 
dimensional integral: 



Gloc(w) = 
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The local self-energy then follows from S]ioc(w) = 
gQ^{oj) — Gi~J(w) which accounts for the implicit (para- 
metric) complex- number relation Eioc[Gioc], i-e. the 
goal is achieved by the semi-analytic exact solution. 
In practice this is done by a parametrization of the 
above integral equation through 5o[Gioc] (inversion), and 
iterating until self-consistency is reached. This works 
fine here only because the interaction lines carry no 
time dependence. In Fig.|2]we show for various disorder 
strengths the local self-energy obtained for E' = 
after convergence, i.e. the answer as predicted by the 
conventional DMFT approach. 

The full calculation involves Monte Carlo sampling of 
all skeleton diagrams except those contributing to Eioc 
(which would otherwise consume about 90% of the sim- 
ulation time already for V = \/2). In the real-space rep- 
resentation this means that only skeleton graphs which 
contain at least two vertices with different site indices are 
accounted for in E'. The simulation itself was done using 
standard BDMC rules with the self-consistency loop im- 
plemented exactly as described in the introductory part 
of this Letter. It turns out that the diagrammatic se- 
ries for Anderson's localization problem constitutes the 
'worst case scenario' in terms of convergence properties. 



Although for any finite time t the series are convergent 
(allowing us to use Dyson's equation and Eq. ([6|), the 
required expansion order increases dramatically with the 
time t. Realistically, we were able to deal with skeleton 
graphs up to order rimax ^ 50 which was limiting the 
accessible times in the simulation of E'. We observe that 
the values of E^_j., = E^ turn out to be extremely small, 
about two orders of magnitude smaller than Ejoc even 
in the intermediate coupling regime V = -^2, see Fig. [s] 
Since the complexity, and hence the relative error-bar, of 
the BDMC simulation is roughly the same for simulating 
E or E', we conclude that the BDMC-hDMFT scheme 
produces results which are two orders of magnitude (or 
a speedup of ^ 10** in CPU time) more accurate for the 
same simulation time in the region of parameter space 
where the series converges and error bars are under con- 
trol. This constitutes the proof of principle for the pro- 
posed scheme. Final results for the density of states are 
indistinguishable from the exact diagonalization data. 




FIG. 2: (Color online) Local self-energy calculated using local 
propagators (S' = 0) for disorder strengths V — l/v^, \/2, 4 
on a lattice of size 12 x 12 x 12. 

Outlook and Conclusions. - We have introduced an 
approach that uses DMFT as an integral part of perform- 
ing simulations of skeleton graphs in strongly interacting 
systems. It combines the power of solving impurity 
problems efficiently with the diagrammatic formalism 
that is unbiased and exact. Given the already good 
agreement between DMFT and diagrammatic Monte 
Carlo based on bare propagators for the Hubbard model 
at U / J — 4 0], we expect the present formalism to bring 
radical speed up and accuracy to studies of the Hub- 
bard model at larger values of U and lower temperatures. 

We would also like to mention several generalizations 
of the simplest scheme introduced above. To begin with, 
the definition of momentum-independent propagator 
allows the use of an arbitrary function /(p) in the defi- 
nition of Gioc such that Gioc = J^^ ^(p)/(p) ^P/(2^)'^- 
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FIG. 3: (Color online) Correction to local self-energy for V = 
a/2 on a lattice of size 12 x 12 x 12: E'jj|^o(t) (red line) and 
^'n|=i(^) (blue line). The noise in the curve is indicative of 
the error bars. The sign problem and the high expansion 
orders put a limit on the accessible times. 



with specific momentum-dependence and structure, (and 
compensating them by impurity solvers dealing with 
a few sites, similar to the ideas behind cluster-DMFT 
schemes). The diagrams to be summed up by the impu- 
rity solver are those with the connections of a compact 
cluster of sites. Similar extensions for real-space clusters 
are also possible. 
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The rest of the scheme remains intact: as before 
diagrams containing exclusively Gioc propagators are 
all summed up in the local self-energy while S' con- 
tains at least one line which is based on G{p) — Gioc- 
The freedom of choosing /(p) different from a constant 
may be used to optimize the subtraction of leading terms. 

In the generic many-body skeleton diagram, any renor- 
malized line whether it is the single-particle propagator 
G{p), the interaction line W^(q), or the two-particle 
propagator F, can be split into momentum-independent 
and momentum-dependent parts (with the same freedom 
of defining the local part as described in the previous 
paragraph). Next, all diagrams based exclusively on 
momentum-independent lines can be dealt with using 
impurity solvers with BDMC accounting for the remain- 
ing graphs. Since the summation of certain geometric 
series such as ladder or screening diagrams can be done 
analytically to set up the original diagrammatic space, 
one can go even further beyond the purely local physics 
by doing so. 

Our final remark is that nothing prevents one 
from extending the idea of subtracting diagrams with 
momentum-independent lines (and compensating them 
separately by impurity solvers) to subtracting diagrams 
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